%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerische Blochgleichungen mit Relaxation, n = (n11 n22 n12 n21).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch

    figure(1);					clf;
    V=[0:.1:2*pi];

    subplot(221),plot3(sin(V),cos(V),0*V,'b','MarkerSize',10,'LineWidth',.01)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)
    axis equal;
    view([-20 20])
    hold on
    subplot(221),plot3(0*V,sin(V),cos(V),'b','MarkerSize',10,'LineWidth',.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)
    subplot(221),plot3(cos(V),0*V,sin(V),'b','MarkerSize',10,'LineWidth',.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)
    subplot(222),plot3(sin(V),cos(V),0*V,'b','MarkerSize',10,'LineWidth',.01)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)
    axis equal;
    view([-20 20])
    hold on
    subplot(222),plot3(0*V,sin(V),cos(V),'b','MarkerSize',10,'LineWidth',.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)
    subplot(222),plot3(cos(V),0*V,sin(V),'b','MarkerSize',10,'LineWidth',.2)
    set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',.01)

    dt  =      3e-9;
    Ga  =  2*pi*0e6;
    De  =  2*pi*6e6;
    Om  =  2*pi*6e6;
    T   =  pi/De*1;

    rho0=[1 0 0 0]';	
    t   =         0;
    while t<T
    	d=0 *(t<T);	o=Om*(t<T);
    	
    	t=t+dt;
        rho = Bloch2(rho0,dt,Ga,o,d);
    	
    	SubPlot(1,t,dt,rho0,rho);
    	rho0=rho;
    end

    rho0=[.5 .5 .5 0]';	
    t   =            0;
    while t<T
        d=De *(t<T);	o=0*(t<T);

        t=t+dt;
        rho = Bloch2(rho0,dt,Ga,o,d);

        SubPlot(2,t,dt,rho0,rho);
        rho0=rho;
    end

    hold off
end


function SubPlot(pic,t,dt,rho0,rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plots.

	n0(1)=2*imag(rho0(3));
	n0(2)=2*real(rho0(3));
	n0(3)=rho0(2)-rho0(1);
	n(1) =2*imag(rho(3));
	n(2) =2*real(rho(3));
	n(3) =rho(2)-rho(1);
	
	subplot(2,2,pic),plot3([n0(1) n(1)],[n0(2) n(2)],[n0(3) n(3)],'.r','MarkerSize',10,'LineWidth',1.8)
	set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',18,'LineWidth',2)
	xlabel('\fontsize{18}Im {\it\rho_{12}}')
	ylabel('\fontsize{18}Re {\it\rho_{12}}')
	zlabel('\fontsize{18}2{\it\rho_{22}}-1')
end


function n = Bloch2(n0,t,g,o,d)
    I =     1i/2;

    M=[ 0      g     I*o   -I*o 
        0     -g	-I*o    I*o 
        I*o  -I*o -g/2-1i*d  0   
       -I*o   I*o     0   -g/2+1i*d];
    n=expm(M*t)*n0;
end
